      subroutine setcfls

      use param

#include <global.h>
      include 'ntbytes.h'
      include 'geom.h'
      include 'xsec.h'
      include 'setls.h'
      include 'sanm.h'
      include 'scmfdp1.h'

! flux/current correction factor
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          do iem=1,2
            dtl(iem,j,im)=2*beta(iem,jm1,im)*beta(iem,j,im)/(beta(iem,jm1,im)+beta(iem,j,im))
            dtr(iem,j,im)=2*beta(iem,j,im)*beta(iem,jp1,im)/(beta(iem,j,im)+beta(iem,jp1,im))
          enddo          
        enddo
        dtl(1,1,im)=bcm_p1(1)
        dtr(1,tnode,im)=bcm_p1(2)
        dtl(2,1,im)=dtl(1,1,im)
        dtr(2,tnode,im)=dtr(1,tnode,im)
      enddo
      
      do im=1,ng
        do j=1,tnode
          jm1=j-1
          jp1=j+1
          btm1=beta(1,jm1,im)
          bt1=beta(1,j,im)
          btp1=beta(1,jp1,im)
          btm2=beta(2,jm1,im)
          bt2=beta(2,j,im)
          btp2=beta(2,jp1,im)
          ppm1=psphi(1,jm1,im)
          pp1=psphi(1,j,im)
          ppp1=psphi(1,jp1,im)
          ppm2=psphi(2,jm1,im)
          pp2=psphi(2,j,im)
          ppp2=psphi(2,jp1,im)

          dhl(1,j,im)=-(jnet(1,j,im)+dtl(1,j,im)*(pp1-ppm1))/(pp1+ppm1)   ! current correction factor
          dhr(1,j,im)=-(jnet(1,jp1,im)+dtr(1,j,im)*(ppp1-pp1))/(ppp1+pp1)
          dhl(2,j,im)=-(jnet(2,j,im)+dtl(2,j,im)*(pp2-ppm2))/(pp2+ppm2)   
          dhr(2,j,im)=-(jnet(2,jp1,im)+dtr(2,j,im)*(ppp2-pp2))/(ppp2+pp2)

          sffdm1=(bt1*pp1+btm1*ppm1)/(bt1+btm1)  ! surface flux for the FDM
          sffdm2=(bt2*pp2+btm2*ppm2)/(bt2+btm2) 
          flcf(1,j,im)=(sflux(1,j,im)-sffdm1)/(pp1+ppm1)  ! flux correction factor
          flcf(2,j,im)=(sflux(2,j,im)-sffdm2)/(pp2+ppm2)
!          btm2=beta(2,jm1,im)
!          bt2=beta(2,j,im)
!          sffdm2=(bt2*avgflx(2,j,im)+btm2*avgflx(2,jm1,im))/(bt2+btm2)
          if (j.eq.tnode) then
            sffdm1=(btp1*ppp1+bt1*pp1)/(btp1+bt1)
            sffdm2=(btp2*ppp2+bt2*pp2)/(btp2+bt2)
            flcf(1,jp1,im)=(sflux(1,jp1,im)-sffdm1)/(pp1+ppp1)
            flcf(2,jp1,im)=(sflux(2,jp1,im)-sffdm2)/(pp2+ppp2)
          endif
        enddo
        if (bcb.eq.1) then
          dhl(1,1,im)=(0.375*sflux(2,1,im)-jnet(1,1,im))/sflux(1,1,im)-dtl(1,1,im)
          dhl(2,1,im)=(0.125*sflux(1,1,im)-jnet(2,1,im))/sflux(2,1,im)-dtl(2,1,im)
        endif
        if (bcu.eq.1) then
          dhr(1,tnode,im)=(-0.375*sflux(2,tnode+1,im)-jnet(1,tnode+1,im))/sflux(1,tnode+1,im) &
                         +dtr(1,tnode,im) 
          dhr(2,tnode,im)=(-0.125*sflux(1,tnode+1,im)-jnet(2,tnode+1,im))/sflux(2,tnode+1,im) &
                         +dtr(2,tnode,im) 
        endif
      enddo

      end subroutine
